version 11.0
set more off

*     ***************************************************************** *;
*     ***************************************************************** *;
*       File-Name:      MiltonPrice_II_Final.do		                	*;
*       Date:           12/18/2019                                      *;
*       Purpose:        To replicate main results in Milton and Price's	*;
*						International Interactions article 				*;
*       Input File:     MiltonPrice_II_Final.dta   		            	*;
*     ****************************************************************  *;
*     ****************************************************************  *;

*     ****************************************************************  *;
*		Table 1, Model 1												*;
*     ****************************************************************  *;

logit failure size religious ethnic democracy civilwar logpop namerica latamerica ssafrica meast asia spline1 spline2 spline3, robust cluster(groupid)

*     ****************************************************************  *;
*		Table 1, Model 2												*;
*     ****************************************************************  *;

logit failure decap size religious ethnic democracy civilwar logpop namerica latamerica ssafrica meast asia spline1 spline2 spline3, robust cluster(groupid)

*     ****************************************************************  *;
*		Table 1, Model 3												*;
*     ****************************************************************  *;

logit failure allies size religious ethnic democracy civilwar logpop  namerica latamerica ssafrica meast asia spline1 spline2 spline3, robust cluster(groupid)

*     ****************************************************************  *;
*		Table 1, Model 4												*;
*     ****************************************************************  *;

logit failure decap allies interact_allies size religious ethnic democracy civilwar logpop  namerica latamerica ssafrica meast asia spline1 spline2 spline3, robust cluster(groupid)

*     ****************************************************************  *;
*		Table 1, Model 5												*;
*     ****************************************************************  *;

logit failure alliesties size religious ethnic democracy civilwar logpop  namerica latamerica ssafrica meast asia spline1 spline2 spline3, robust cluster(groupid)

*     ****************************************************************  *;
*		Table 1, Model 6												*;
*     ****************************************************************  *;

logit failure decap alliesties interact size religious ethnic democracy civilwar logpop namerica latamerica ssafrica meast asia spline1 spline2 spline3, robust cluster(groupid)

*     ****************************************************************  *;
*     ****************************************************************  *;
*     ****************************************************************  *;
*          Run This Code Through to End to Create Figure 1              *;
*     ****************************************************************  *;
*     ****************************************************************  *;
*     ****************************************************************  *;           
        
#delimit ;
		
logit failure decap alliesties interact size religious ethnic democracy civilwar logpop namerica latamerica ssafrica meast asia spline1 spline2 spline3, robust cluster(groupid);
        
*     ****************************************************************  *;
*       Take 10,000 draws from the estimated coefficient vector and     *;
*       variance-covariance matrix.                                     *;
*     ****************************************************************  *;

preserve;

set seed 339487731;

drawnorm SN_b1-SN_b18, n(10000) means(e(b)) cov(e(V)) clear;

*     ****************************************************************  *;
*       To calculate the desired quantities of interest we need to set  *;
*       up a loop.  This is what we do here.                            *;
*     ****************************************************************  *;
*       First, specify what you quantities should be saved and what     *;
*       these quantities should be called.                              *;
*     ****************************************************************  *;

postutil clear;
postfile mypost prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi 
            using sim.dta, replace;
            
noisily display "start";
            
*     ****************************************************************  *;
*       Start loop.  Let `a' be the modifying variable Z and let this   *;
*       run from min to max in the desired increments.                  *;
*     ****************************************************************  *; 

local a=0 ;
while `a' <= 1 { ;


{;

scalar h_decap=0;
scalar h_alliesties=0.0325836;
scalar h_size=1;
scalar h_religious=1;
scalar h_ethnic=0;
scalar h_democracy=1;
scalar h_civilwar=0;
scalar h_logpop=15.86592;
scalar h_namerica=0;
scalar h_latamerica=0;
scalar h_ssafrica=0;
scalar h_meast=1;
scalar h_asia=0;
scalar h_spline1=-340.8906;
scalar h_spline2=-444.148;
scalar h_spline3=-338.1351;
scalar h_constant=1;          

    generate x_betahat0  = SN_b1*h_decap 
                            + SN_b2*`a'   
                            + SN_b3*h_decap*`a' 
                            + SN_b4*h_size 
                            + SN_b5*h_religious
                            + SN_b6*h_ethnic
                            + SN_b7*h_democracy
                            + SN_b8*h_civilwar
							+ SN_b9*h_logpop
							+ SN_b10*h_namerica
							+ SN_b11*h_latamerica
							+ SN_b12*h_ssafrica
							+ SN_b13*h_meast
							+ SN_b14*h_asia
							+ SN_b15*h_spline1
							+ SN_b16*h_spline2
							+ SN_b17*h_spline3
                            + SN_b18*h_constant;
                            
    generate x_betahat1  = SN_b1*(h_decap+1) 
                            + SN_b2*`a'   
                            + SN_b3*(h_decap+1)*(`a') 
                            + SN_b4*h_size 
                            + SN_b5*h_religious
                            + SN_b6*h_ethnic
                            + SN_b7*h_democracy
                            + SN_b8*h_civilwar
							+ SN_b9*h_logpop
							+ SN_b10*h_namerica
							+ SN_b11*h_latamerica
							+ SN_b12*h_ssafrica
							+ SN_b13*h_meast
							+ SN_b14*h_asia
							+ SN_b15*h_spline1
							+ SN_b16*h_spline2
							+ SN_b17*h_spline3
                            + SN_b18*h_constant; 

    gen odds0=exp(x_betahat0);
    gen odds1=exp(x_betahat1);
							
    gen prob0=odds0/(1+odds0);
    gen prob1=odds1/(1+odds1);
    gen diff=prob1-prob0;
    
    egen probhat0 =mean(prob0);
    egen probhat1=mean(prob1);
    egen diffhat=mean(diff);
    
    tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi;  

    _pctile prob0, p(5,95) ;
    scalar `lo0' = r(r1);
    scalar `hi0' = r(r2);  
    
    _pctile prob1, p(5,95);
    scalar `lo1'= r(r1);
    scalar `hi1'= r(r2);  
    
    _pctile diff, p(5,95);
    scalar `diff_lo'= r(r1);
    scalar `diff_hi'= r(r2);  
   
    scalar `prob_hat0'=probhat0;
    scalar `prob_hat1'=probhat1;
    scalar `diff_hat'=diffhat;
    
    post mypost (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') (`lo1') (`hi1') 
                (`diff_hat') (`diff_lo') (`diff_hi');
   
    };      
    
    drop x_betahat0 x_betahat1 odds0 odds1 prob0 prob1 diff probhat0 probhat1 diffhat;
    
    local a=`a'+ .01;
    
    display "." _c;
    
} ;

display "";

postclose mypost;

*     ****************************************************************  *;                                  
*       Call on posted quantities of interest                           *;
*     ****************************************************************  *;

restore;

merge using sim.dta;

gen yline=0;

gen MV = (_n-1)/100;

replace  MV=. if _n>65;

#delimit ;
			
graph twoway hist alliesties if alliesties<.16, percent color(gs14) yaxis(2)
        ||  line diff_hat MV if MV<.16, clpattern(solid) clwidth(medium) clcolor(blue) clcolor(black)
        ||  line diff_lo MV if MV<.16, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line diff_hi MV if MV<.16, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line yline MV if MV<.16,  clwidth(thin) clcolor(black) clpattern(solid)
        ||  ,   
            xlabel(, nogrid labsize(2)) 
            ylabel(, axis(1) nogrid labsize(2))
            ylabel(, axis(2) nogrid labsize(2))
            yscale(noline alt) yscale(noline alt axis(2)) xscale(noline) legend(off)
            yline(0, lcolor(black))
            xtitle("Allies' Ties", size(2.5))
            ytitle("Marginal Effect of Leadership Decapitation", size(2.5))
            ytitle("Percent of Observations" , axis(2) size(2.5))
            xsca(titlegap(4)) ysca(titlegap(4)) ysca(axis(2) titlegap(4))
            scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));

stop; 			

*Running this code will plot the marginal effect across all values of allies' ties		

#delimit ;
			
graph twoway hist alliesties, percent color(gs14) yaxis(2)
        ||  line diff_hat MV, clpattern(solid) clwidth(medium) clcolor(blue) clcolor(black)
        ||  line diff_lo MV, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line diff_hi MV, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line yline MV,  clwidth(thin) clcolor(black) clpattern(solid)
        ||  ,   
            xlabel(, nogrid labsize(2)) 
            ylabel(, axis(1) nogrid labsize(2))
            ylabel(, axis(2) nogrid labsize(2))
            yscale(noline alt) yscale(noline alt axis(2)) xscale(noline) legend(off)
            yline(0, lcolor(black))
            xtitle("Allies' Ties", size(2.5))
            ytitle("Marginal Effect of Leadership Decapitation", size(2.5))
            ytitle("Percent of Observations" , axis(2) size(2.5))
            xsca(titlegap(4)) ysca(titlegap(4)) ysca(axis(2) titlegap(4))
            scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));

			


